/* Copyright 2007 dnAnalytics Project.
 *
 * Contributors to this file:
 * Marcus Cuda
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

#ifndef LINEAR_ALGEBRA_H	


#ifdef  __cplusplus
	extern "C" {
	#endif	
	
	size_t s_index_max(const int n, const float data[]);
	size_t d_index_max(const int n, const double data[]);	
	size_t c_index_max(const int n, const MKL_Complex8 data[]);
	size_t z_index_max(const int n, const MKL_Complex16 data[]);

	float s_absolute_sum(const int n, const float data[]);
	double d_absolute_sum(const int n, const double data[]);
	float c_absolute_sum(const int n, const MKL_Complex8 data[]);
	double z_absolute_sum(const int n, const MKL_Complex16 data[]);

	float s_dot_product(const int n, const float x[], const float y[]);
	double d_dot_product(const int n, const double x[], const double y[]);
	MKL_Complex8 c_dot_product(const int n, const MKL_Complex8 x[], const MKL_Complex8 y[]);
	MKL_Complex16 z_dot_product(const int n, const MKL_Complex16 x[], const MKL_Complex16 y[]);
	MKL_Complex8 c_dot_product_c(const int n, const MKL_Complex8 x[], const MKL_Complex8 y[]);
	MKL_Complex16 z_dot_product_c(const int n, const MKL_Complex16 x[], const MKL_Complex16 y[]);

	void s_scale(const int n, const float alpha, float x[]);
	void d_scale(const int n, const double alpha, double x[]);
	void c_scale(const int n, const MKL_Complex8 alpha, MKL_Complex8 x[]);
	void z_scale(const int n, const MKL_Complex16 alpha, MKL_Complex16 x[]);

	void s_axpy( const int n, const float alpha, const float x[], float y[]);
	void d_axpy( const int n, const double alpha, const double x[], double y[]);
	void c_axpy( const int n, const MKL_Complex8 alpha, const MKL_Complex8 x[], MKL_Complex8 y[]);
	void z_axpy( const int n, const MKL_Complex16 alpha, const MKL_Complex16 x[], MKL_Complex16 y[]);

	void d_vector_add( const int n, const double x[], const double y[], double ret[]);
	void d_vector_subtract( const int n, const double x[], const double y[], double ret[]);
	void d_vector_multiply( const int n, const double x[], const double y[], double ret[]);
	
	void s_matrix_multiply(const int m, const int n, const int k, const float x[],
		const float y[], float c[]);
	void d_matrix_multiply(const int m, const int n, const int k, const double x[],
		const double y[], double c[]);
	void c_matrix_multiply(const int m, const int n, const int k, const MKL_Complex8 x[],
		const MKL_Complex8 y[], MKL_Complex8 c[]);
	void z_matrix_multiply(const int m, const int n, const int k, const MKL_Complex16 x[],
		const MKL_Complex16 y[], MKL_Complex16 c[]);

	int s_lu_factor(int m, float a[], int ipiv[]);
	int d_lu_factor(const int m, double a[], int ipiv[]);
	int c_lu_factor(int m, MKL_Complex8 a[], int ipiv[]);
	int z_lu_factor(int m, MKL_Complex16 a[], int ipiv[]);

	int s_lu_inverse(int n, float a[], int ipiv[]);
	int d_lu_inverse(int n, double a[], int ipiv[]);
	int c_lu_inverse(int n, MKL_Complex8 a[], int ipiv[]);
	int z_lu_inverse(int n, MKL_Complex16 a[], int ipiv[]);

	int s_lu_solve(int n, int nrhs, float a[], int ipiv[], float b[]);
	int d_lu_solve(int n, int nrhs, double a[], int ipiv[], double b[]);
	int c_lu_solve(int n, int nrhs, MKL_Complex8 a[], int ipiv[], MKL_Complex8 b[]);
	int z_lu_solve(int n, int nrhs, MKL_Complex16 a[], int ipiv[], MKL_Complex16 b[]);

	float s_matrix_norm(char norm, int m, int n, float a[]);
	double d_matrix_norm(char norm, int m, int n, double a[]);
	float c_matrix_norm(char norm, int m, int n, MKL_Complex8 a[]);
	double z_matrix_norm(char norm, int m, int n, MKL_Complex16 a[]);

	int s_condition(char type, int n, float norm, float a[], float* rcond);
	int d_condition(char type, int n, double norm, double a[], double* rcond);
	int c_condition(char type, int n, float norm, MKL_Complex8 a[], float* rcond);
	int z_condition(char type, int n, double norm, MKL_Complex16 a[], double* rcond);

	int s_cholesky_factor(int n, float a[]);
	int d_cholesky_factor(int n, double a[]);
	int c_cholesky_factor(int n, MKL_Complex8 a[]);
	int z_cholesky_factor(int n, MKL_Complex16 a[]);

	int s_cholesky_solve(int n, int nrhs, float a[], float b[]);
	int d_cholesky_solve(int n, int nrhs, double a[], double b[]);
	int c_cholesky_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[]);
	int z_cholesky_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[]);

	int s_qr_factor(int m, int n, float r[], float tau[], float q[]);
	int d_qr_factor(int m, int n, double r[], double tau[], double q[]);
	int c_qr_factor(int m, int n, MKL_Complex8 r[], MKL_Complex8 tau[], MKL_Complex8 q[]);
	int z_qr_factor(int m, int n, MKL_Complex16 r[], MKL_Complex16 tau[], MKL_Complex16 q[]);

	int s_qr_solve(int m, int n, int bn, float r[], float b[], float tau[], float x[]);
	int d_qr_solve(int m, int n, int bn, double r[], double b[], double tau[], double x[]);
	int c_qr_solve(int m, int n, int bn, MKL_Complex8 r[], MKL_Complex8 b[], MKL_Complex8 tau[], MKL_Complex8 x[]);
	int z_qr_solve(int m, int n, int bn, MKL_Complex16 r[], MKL_Complex16 b[], MKL_Complex16 tau[], MKL_Complex16 x[]);
	
	int s_svd_factor(bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[]);
	int d_svd_factor(bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[]);
	int c_svd_factor(bool compute_vectors, int m, int n, MKL_Complex8 a[], float s[], MKL_Complex8 u[], MKL_Complex8 v[]);
	int z_svd_factor(bool compute_vectors, int m, int n, MKL_Complex16 a[], double s[], MKL_Complex16 u[], MKL_Complex16 v[]);
	
	double d_householder(int n, double x[]);

	void d_transpose(int row, int col, double a[], double b[]);

	#ifdef  __cplusplus
		}
	#endif

	#define LINEAR_ALGEBRA_H

#endif